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. , ABSTRACT 

■ We use a model for the evolution of galaxies in the far-IR based on the ACDM cosmology to 

make detailed predictions for upcoming cosmological surveys with the Herschel Space Obser- 

Q . vatory. We use the combined GALFORM semi-analytical galaxy formation model and GRAS I L 

|H ' spectrophotometric code to compute galaxy SE Ds including the repr ocessing of radiation by 

C/5 \ dust. The model, which is the same as that in iBaugh et al.l (l2005b . assumes two different 

03 . IMFs: a normal solar neighbourhood IMF for quiescent star formation in disks, and a very 

top-heavy IMF in starbursts triggered by galaxy mergers. We have shown previously that the 

, top-heavy IMF appears necessary to explain the number counts and redshifts of faint sub-mm 

^ ' galaxies. In this paper, we present predictions for galaxy luminosity functions, number counts 

I and redshift distributions in the Herschel imaging bands. We find that source confusion will 

\^ , be a serious problem in the deepest planned surveys. We also show predictions for physical 

•/^ ' properties such as star formation rates and stellar, gas and halo masses, together with fluxes 

\ at other wavelengths (from the far-UV to the radio) relevant for multi-wavelength follow-up 

Q>^ i observations. We investigate what fraction of the total IR emission from dust and of the high- 

' mass star formation over the history of the Universe should be resolved by planned surveys 

, with Herschel, and find a fraction ^ 30 — 50%, depending on confusion. Finally, we show 

■ that galaxies in Herschel surveys should be significantly clustered. 
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1 INTRODUCTION 

The Herschel Space Observatory was launched on 14 May 2009 
and will begin science observations towards the end of the year. It 
will observe the Universe at far-infrared (IR) wavelengths, from 60 
to 670/im, and will be far more sensitive at these wavelengths than 
any previous telescope. One of its primary goals will be to probe the 
dust-obscured part of the cosmic history of star formation. In this 
paper we use a state-of-the-art theoretical model of galaxy forma- 
tion based on structure formation in the ACDM model, combined 
with a detailed treatment of the reprocessing of stellar radiation to 
far-IR wavelengths, to make predictions for what should be seen 
in cosmological surveys with Herschel, and for how well Herschel 
should be able to achieve its goal of unveiling the cosmic star for- 
mation history. 

The Herschel satellite follows in a line of previous space- 
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based IR telescopes, starting in the 1980s and 1990s with IRAS 
and ISO and including most recently Spitzer and AKARI. These 
have gradually revealed the nature and evolution of the galaxy 
population at mid- and far-IR wavelengths which are inaccessi- 
ble from the ground, where galaxy luminosities are d ominated by 
emission from interstellar dust grains (see reviews bv lSoifer et al.l 
1 19871 : lElbazll2005l : ISoifer et al.ll200Sl . More recently, these obser- 
vations from space have been complemented by surveys at longer, 
sub-mm wavelengths from the ground, starting wit h surveys at 
850^m using the SCUBA instrument on the JCMT jSmail et"ail 
ll997l : lHughes et al.lll998l) . However, due to their poor sensitivities 
and angular resolutions at these wavelengths, these earlier space 
missions provided only very limited direct information on the evo- 
lution of galaxies in the rest-frame far-IR range which contains 
most of the energy re-radiated from dust grains heated by starlight. 
The ground-based sub-mm surveys have also been limited to seeing 
only the very highest luminosity galaxies at high redshifts, and have 
been hampered by the fact that they only observe dust emission 
longwards of the peak in the spectral energy distribution (SED). 
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Herschel will, for the first time, allow observations of galaxy SEDs 
around the far-IR peak in dust emission out to high redshifts, and 
thus back to when the Universe was only a fraction of its current 
age. 

Although observations with IRAS had already shown that cer- 
tain types of nearby star-forming galaxies (the Ultra-Luminous IR 
Galaxies, or ULIRGs) emit most of thei r luminosity thro ugh dust 
emission in the far-IR (see the review bY Soiferetal.ll987h . a land- 
mark was achie ved with the discovery by COBE of the cosmic far- 
IR ba ckground jPuget et al1ll996l : lHa"user et al.|[T993 : lFixsen et alj 
Il998l) . This far-IR background was found to have an energy density 
comparable to the ultraviolet (UV)/optical background from stars. 
The far-IR background is most naturally interpreted as the emis- 
sion from dust in galaxies heated by starlight, integrated over the 
history of the Universe. (There is also a contribution to the far-IR 
background from dust heated by radiation from AGN, but, based on 
multi-wavelength, especially X-ray, studies, t he AGN contribution 
appears to be sma ll, only ^10% overall (e.g. lAlmaini et al. I I 19991 : 
iFardal et alj|2007t) . Since most of the heating of dust in galaxies 
is due to radiation from young stars, the far-IR background pro- 
vides powerful evidence that the bulk of star formation over the 
history of the Universe has been obscured by dust, with most of 
the radiation from young stars having been reprocessed from UV 
to far-IR wavelengths. It is hoped that observations with Herschel 
will allow most of the far-IR background to be directly resolved 
into individual sources. Preliminary constraints on the far-IR source 
counts ha ve recently been ob tained with the balloon-born BLAST 
telescope jPevUn et aLll2009h . but surveys with Herschel will have 
both better angular resolution and much better sensitivity. 

Measuring the cosmic star formation history, and understand- 
ing it in terms of physical models for galaxy formation and evo- 
lution, are among the main goals of modem cosmology. Obser- 
vational studies of the cosmic SFR history began with optically- 
selected samples of galaxies at different redshifts, which typically 
used th e rest-frame UV emission from galaxi es as the SFR in - 
dicator juillv et al.ltl994 iMadau et al.l 119961 : Isteidel et al.lll999l) . 
The earliest studies derived cosmic SFR histories ignoring the ef- 
fects of dust extinction altogether, but it soon became apparent that 
this approach was inadequate, since studies of local star-forming 
galaxies had previously shown that these galaxies us ually had ap- 
preci able dust extinctions in the rest-frame UV (e.g. iMeurer et al.l 
1 19951) . Applying locall y-derived relations between UV extinct ion 
and UV spectral slope jlVIeurer et al]| 19951 : Icalzetti et alJI 19951) to 
high-redshift optically-selected galaxy samples, in particular the 
Lyman-break galaxies (LBGs), implied that these galaxies should 
have large UV dust extinctions, and thus also large dust correc- 
tions to the SFRs i nferred from their rest-frame UV luminosities 
jSteidel et al.ll999l : lMeurer et aljl999l) . However, in the absence of 
rest-frame far-IR data on these high-redshift galaxies, which would 
directly measure the amount of energy absorbed and re-radiated 
by dust, these corrections for dust extinction remain very uncer- 
tain. Correspondingly the cosmic SFR histories derived from dust- 
corrected rest-frame UV data also have large uncertainties. 

There have also been atte mpts to estimate the cosmic SFR his- 
tory f rom sub-mm surveys (e.g. lHughes et al.ll998l : IChaDman et al.l 
l2005h . but here the problems have been that these surveys have de- 
tected only small numbers of the most luminous galaxies at high 
redshift, that the redshifts of many of these are uncertain, and that 
an extrapolation of the SED must be made from the sub-mm to 
the far-IR in order to derive the total IR luminosity from dust, 
from which the SFR is calculated. (In a highly dust-obscured star- 
forming galaxy, almost all of the UV light from young stars is re- 



processed into IR emission by dust, and so the total dust luminosity 
provides a good measure of the underlying UV luminosity, and thus 
of the SFR.) More recently, observations of the mid-IR emission 
from dust in galaxies have been used to try to infer the cosmic SFR 
history (e.g. lUe Floc'h et al.l200^ : |Perez-Gonzalez et alj2005h , but 
this method has the drawback that the SED must be extrapolated 
from the mid-IR to the far-IR in order to estimate the total dust 
luminosity. The mid-IR/far-IR ratio is known from nearby exam- 
ples to show large variations between galaxies, so this extrapolation 
is likewise uncertain when applied to high-redshift galaxies where 
there is no direct measurement of the far-IR. Measurements of the 
cosmic SFR history using Herschel will avoid all of the problems 
associated with these other SFR tracers by measuring the far-IR 
emission directly. 

An important issue which we have not yet discussed is that 
of the stellar Initial Mass Function (IMF). All of the methods 
used to estimate observationally the SFRs of high-redshift galaxies 
(whether from rest-frame UV or IR emission, or emission line or 
radio luminosities) actually only directly trace massive young stars 
(typically > 5Mq), and so really provide measures of the SFR for 
high-mass stars only. Deriving the total SFR for the whole range 
of stellar mass (~ 0.1 — lOOi\f0) requires an extrapolation using 
an assumed IMF. It is conventional in observational estimates of 
the cosmic SFR history to assume a universal IMF, usually taken 
to be similar to that in the solar neighbourhood. However, this as- 
sumption of a universal IMF has no direct observational basis, es- 
pecially at higher redshifts, and indeed already appears to lead to 
inconsistencies between estimated cosmic SFR histories and inde- 
pendent estimates of the evolution of the stellar mass density (e.g. 
Perez-Gonzalez et al. 2008). We return to this issue later in the pa- 
per. 

As already stated, the aim of this paper is to make predic- 
tions for what cosmological surveys with Herschel should see, 
based on the best current ab initio models both for how galax- 
ies form and evolve, and for how the radiation they emit is dis- 
tributed over UV, optical, IR and sub-mm wavelengths. Before 
we discuss our own approach in more detail, we briefly review 
the different modelling strategies for galaxy evolution in the mid- 
and far-IR. The models can be broadly divided into two classes, 
phenomenological models and hierarchical galaxy formation mod- 
els. In phenomenological models, the galaxy luminosity function 
(LF) in the IR and its evolution are given by some phenomeno- 
logical fit, and the IR SEDs are likewise empirical fits, with the 
parameters in the evolving LF adjusted to match various obser- 
vational data (e.g. Rowan-Robinson 200 ll: iLagache et al] l2003l : 
lRowan-Robinsonll2009l) . These models are generally restricted to 
modelling the dust emission from galaxies, and do not include the 
UV /optical/near-IR emission from stars. They have been used by 
e.g. iFemandez-Conde et al.l ( 1200 Sl) to make predictions for Her- 
schel. In hierarchical galaxy formation models, the mass assembly 
of galaxies is related to structure formation in the dark matter, and 
star formation and galaxy merger histories are calculated based on 
physical prescriptions for star formation, supernova feedback, dy- 
namical friction etc. These models typically use theoretical stellar 
population synthesis models to compute galaxy stellar luminosi- 
ties, combined with some physical model for the dust extinction, 
but then can be distinguished according to how they compute the 
SED of dust emission. Som e employ empirical or phenom elogi- 
cal SEDs for the dust (e.g. lOevriendt & Guiderdonil |2000|) : oth- 
ers employ fully theoretical dust SEDs based on radiative trans- 
fer and detailed modelling of heating and cooling of dust grains. 
The modelling approach we use here, which has been set out in 
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iGranato et al] feOOOh . ISaugh et alj feOOSi) and lLacev et d] j2008l) . 
is of the latter type, with detailed physical modelling both of galaxy 
formatio n and of the galaxy S EDs. Other models of the same type 
include Fontanot et alj ( |2007|) . We note here some of the advan- 
tages of hierarchical galaxy formation models over phenomenolog- 
ical models for galaxy evolution in the IR: they provide a direct link 
to theoretical models of structure formation; typically they predict 
galaxy SEDs over a wider wavelength range than the IR, includ- 
ing at least the UV and optical; generally they predict a wide range 
of other galaxy properties, such as masses, sizes, gas fractions and 
morphologies; and finally, they allow direct predictions of galaxy 
clustering, without any further assumptions. 

This present paper is the fourth in a series, where we com- 
bine the GALFO RM semi-analytical model of galaxy formation 
( IColeetalj|2000l) with the GRASIL model for stellar and dust 
emission from galaxies jSilva et alll 19981) . The GALFORM model 
computes the evolution of galaxies in the framework of the ACDM 
model of structure formation, based on physical treatments of the 
assembly of dark matter halos by merging, gas cooling in halos, 
star formation and supernova feedback, galaxy mergers, and chem- 
ical enrichment. The GRAS IL model computes the SED of a model 
galaxy from the far-UV to the radio, based on theoretical models of 
stellar evolution and stellar atmospheres, radiative transfer through 
a two-phase dust medium to calculate both the dust extinction and 
dust heating, and a distribution of dust temperatures in each galaxy 
calculated from a detailed dust grai n model. 

In the first paper in the series iGranato et al.ll200(]h . we mod- 
elled the IR properties of galaxies in the local Universe. While this 
model was very successful in explaining observations of the lo- 
cal Universe, we found subsequently that it failed when confronted 
with observations of star-forming galaxies at high redshifts, pre- 
dicting far too few sub-mm galaxies (SMGs) at 2: ~ 2 and Lyman- 
break galaxies (LB Gs) at 2 ~ 3. Therefore, in the second paper 
jBaugh et alj2005l) . we proposed a new version of the model which 
assumes a top-heavy IMF in starbursts (with slope a; = 0, com- 
pared to Salpeter slope x — 1.35), but a normal solar neighbour- 
hood IMF for quiescent star formation. In this new model, the star 
formation parameters were also changed to force more star for- 
mation to happen in bursts. This revised model agreed well with 
both the number counts and redshift distributions of SMGs detected 
at 850/im, and with the rest-frame far-UV luminosity function of 
LBGs at z; ~ 3, while still maintaining consistency with galaxy 
properties in the local Universe such as the optical, near-IR and 
far-IR luminosity functions, and gas fractions, metallicities, mor- 
phologies and sizes. In the third paper ^Lacey et al. 2008) . we com- 
pared predictions from the same iBaugh et alj feOOSh model (with- 
out changing any parameters) with observational data on galaxy 
evolution in the mid-IR from Spitzer and found generally good 
agreement. 

This same model was found bv lLe Delhou et al.l ( |2005|.|2006|) 
and lOrsi et al.l ( l2008h to provide a good match to the observed evo- 
lution of the population of Lya-emitting galaxies, and of their clus- 
tering, over the redshift range 2: ~ 3 — 6. Support for the con- 
troversial assumption of a top-heavy IMF in burs ts came from the 
studies of chemical enrichment in this model by iNagashima et al.l 
who found that the metallicities of both the intergalac- 
tic gas in galaxy clusters and the stars in elliptical galaxies were 
predicted to be significantly lower than observed values if a normal 
IMF was assumed for all star formation, but agreed much better if 
a top-heavy IMF in bursts was assumed, as in .Baugh et aU 

We also mention here an alternative semi-analytical approach 
using the GRASIL SED model which has been developed by 



IGranato et all ^2004) and lSilva et aP j2005l) . This approach differs 
from that in the present paper, in that the treatment of mass as- 
sembly of galaxies is greatly simplified, neglecting halo and galaxy 
mergers, and modelling the disk population phenomenologically, 
but it does include a detailed treatment of the relation between QSO 
and spheroid formation, including feedback from the QSO phase on 
galaxy formation . Predictions for Her sche! from this model have 
been presented in iNegrelloetal] ( 12007) . 

The outline of this paper is as follows: In ij2l we review the 
physical ingredients in our model and the motivations for some of 
the key parameter choices. In [J3]and ^ we present model pre- 
dictions for the evolution of the galaxy luminosity function in the 
Herschel bands, and for the consequent number counts and redshift 
distributions. Next, in !j5] we show predictions for some other key 
physical properties of //er.sc/;e/-selected galaxies, while in [j6l we 
show predictions for fluxes at other wavelengths from the far-UV 
to the radio. In fj7]we try to answer the key cosmological question: 
what fraction of the dust-obscured star formation over the history 
of the Universe should the planned surveys with Herschel be able 
to uncover? In ijs] we briefly discuss what the model predicts for 
clustering of galaxies in the Herschel bands. Finally, we present our 
conclusions in i!9| 



2 MODEL 

The model used for the predicti ons in this paper is iden tical to that 
described in iBaugh et alj j2005l) and lLacevet^ ( I2OO8I) (apart from 
minor code updates), so we give only a very brief summary here. 
We use the GALFORM semi-analytical galaxy formation model to 
predict the physical properties of the galaxy population at differ- 
ent redshifts, and combine it with the GRASIL spectrophotometric 
model to predict the detailed SEDs of model galaxies (including 
emission from dust). 



2.1 GALFORM galaxy formation model 

We compute the formation and evolution of galaxies within the 
framework of the ACDM model of structure formation using the 
semi-analytical galaxy formation model GALFORM. The general 
methodology and approximations b ehind the GALFORM m odel are 
set out in detail in lCole et alj j2000l) (see also the review bv lBaughl 
2006). In summary, the GALFORM model follows the main pro- 
cesses which shape the formation and evolution of galaxies. These 
include: (i) the collapse and merging of dark matter halos; (ii) the 
shock-heating and radiative cooling of gas inside dark halos, lead- 
ing to the formation of galaxy disks; (iii) quiescent star formation 
in galaxy disks; (iv) feedback both from supernova explosions and 
from photo-ionization of the IGM; (v) chemical enrichment of the 
stars and gas; (vi) galaxy mergers driven by dynamical friction 
within common dark matter halos, leading to the formation of stel- 
lar spheroids, and also triggering bursts of star formation. The end 
product of the calculations is a prediction of the numbers and prop- 
erties of galaxies that reside within dark matter haloes of different 
masses. The model predicts the stellar and cold gas masses of the 
galaxies, along with their star formation and merger histories, their 
disk and bulge sizes, and their metallicities. 

The prescriptions and parameters for the different processes 
which we use in this paper are identica l to those adopted by 
iBaugh et al.l j2005l) and iLacev et ai] ( I2OO8I) . The background cos- 
mology is a spatially flat CDM universe with a cosmological con- 
stant, with "concordance" parameters fim = 0.3, f^A = 0.7, 
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fib = 0.04, mdh = i/o/ClOOkms^^Mpc"^) = 0.7. The am- 
plitude of the initial spectram of density fluctuations is set by the 
r.m.s. linear fluctuation in a sphere of radius 8/i~^Mpc, as = 0.93. 
The assembly and merger histories of dark matter halos are com- 
puted using a Mo nte Carlo method based on the extended Press- 
Schechter model dCole et al.ll2OO0l) . When used in galaxy forma- 
tion models, these Monte Carlo merger histories have been found 
to give very similar results to using ha lo merger historie s extracted 
directly from N-body simulations (e.g. lHellv et alj2003h . The evo- 
lution of the baryons within the dark halos is then calculated using 
analytical prescriptions for gas cooling, star formation, supernova 
feedback etc. The parameters of the GALFORM model describing 
these baryonic processes were chosen to reproduce a range of prop- 
erties of present-day galaxies (optical, near-IR and far-IR luminos- 
ity functions, fraction of spheroid- vs disk-dominated galaxies, and 
galaxy disk sizes, gas fractions and metallicities as a function of 
luminosity), as well as the number counts and redshift distribution 
of sub-mm galaxies, and the rest-frame fa r-UV luminosity fu nction 
of Lyman-break galaxies at high redshift teaugh et alj|2005h . 

An important feature of our model is the existence of two 
modes of star formation, the "quiescent" mode in galaxy disks, and 
the "burst" mode triggered by galaxy mergers. In the current model, 
starbursts are triggered by both major galaxy mergers (which trans- 
form stellar disks to spheroids) and minor galaxy mergers (which 
leave the stellar disk of the larger galaxy unchanged). The two 
modes of star formation are assumed to have different stellar Ini- 
tial Mass Functions (IMFs). Both IMFs are taken to be piece- 
wise power laws, with slopes x defined by dN/dlnm oc m~^, 
where A'' is the number of stars and m is the stellar mass (so the 
Salpeter slope is a; = 1.35), and covering a stellar mass range 
0.15 < m < 12OM0. Quiescent star formation in galaxy disks 
is a ssumed to have a solar neighbourhood IMF, for which we use 

the iKermicutl ( Il983h parametrization, wit h slope x = 0.4 for 

m < Mq and a; = 1.5 for m > Mq. (The lKennicut^ ( Il983h IMF 
is similar to other popular parametrizations of the solar neighbour- 
hood IMF, such as that of lKroupai I2OOII) .) Bursts of star formation 
triggered by galaxy mergers are assumed to form stars with a top- 
heavy IMF with slope x = 0. The different IMFs result in different 
luminosities and SEDs for a stellar population, as well as differ- 
ent overall rates of gas return and metal ejection from dying stars. 
These effects are all taken into account self-consistently, based on 
the predictions from stellar evolution models. 

As discussed in detail in iBaugh et al.l ( l2005h , the top-heavy 
IMF in bursts was found to be required in order to reproduce the 
observed number counts and redshift distributions of the faint sub- 
mm galaxies. The top-heavy IMF results both in higher bolometric 
luminosities for young stellar populations, and greater production 
of heavy elements and hence also dust, both effects being important 
for reproducin g the properties o f SMGs in the model. Furthermore, 
as shown by iNagashima et al.l ( l2005al fbl). the predicted chemical 
abundances of the X-ray emitting gas in galaxy clusters and of the 
stars in elliptical galaxies also agree better with observational data 
in a model with the top-heavy IMF in bursts, rather than a universal 
solar neighbourhood IMF. Subsequent work using the same model 
has also shown that it predicts galaxy evoluti on in the mid-IR i n 
good agreement with observations by Spitzer jLacev et alj|2008h . 
A more detailed comparison of the model with the properties of 
observed SMGs has been carried out by ,Swi nbank et al . (2008). As 
shown bv lLe Defliou et alj ( |2005| , |2006|) and lOrsi et all j2008i) . the 
same model also reproduces the observed evolution of the luminos- 
ity func tion and clustering of L yg emitting galaxies at high redshift. 
Finally, lOonzalez et alj j2009l) have made a detailed comparison of 



the model with the luminosity function, colours, morphologies and 
sizes of galaxies in the SDSS survey of the local Universe. 

A variety of other observational evidence has accumulated 
which suggests that the IMF in some environments may be top- 
heavy compared to the solar neighbourhood IM F (see lElmegreenI 
( l2009h for a recent review). iRieke et al] ( Il993h argued for a top- 
heavy IMF in the nearby s tarburst M82, bas ed on modelling its 
integrated properties, while |PM"ra et al.l | |200 7) found possible evi- 
dence for a top-heavy IMF in the ultra-luminous starburst Arp220 
from the relative numbers of supernovae of different types observed 
at radio wavelengths. Evidence has been found for a top-heavy IMF 
in s ome star clusters in int ensely star-forming regions, both in M82 
(e.g.lMcCradv et al.'2003'). and i n our own Galaxy (e.g. lFiger et al.l 
I1999I : IStolteetal.. .2005. ; Haravama et al.|[200a) . Observations of 
both the old and young stellar populations in the central 1 pc of 
our Galaxy also favour a top-heavy IMF (Paumard et al." 200^; 
iManess et al. 2007). In the local Universe,^Meurer et al. ( 2009) find 
evidence for variations in the IM F between galaxies from variations 
in the Ha/UV luminosity ratio. iFardal et al.l ( l2007h found that rec- 
onciling measurements of the optical and IR extragalactic back- 
ground with measurements of the cosmic star formation history 
also s eemed to require an average IMF that was somewhat top- 
heavy. IPerez-Gonzalez et al.l ( I2OO8I) compared observational con- 
straints on the evolution of the star formation rate density and stel- 
lar mass density over cosmic time, and found that reconciling these 
two types of data also favours a more top-heavy IMF at higher red- 
shifts, as had been hinted at by earlier studies. Finally, v an DokkumI 
COOSi) found that reconciling the colour and luminosity evolution 
of early-type galaxies in clusters also favoured a top-heavy IMF. 
iLarsorj ( Il998h summarized other evidence for a top-heavy IMF 
during the earlier phases of galaxy evolution, and argued that this 
could be a natural consequence of the temperature-depend ence of 
the Je ans mass for gravitational instability in gas clouds. iLarsonl 
( I2OO5I) extended this to argue that a top-heavy IMF might also be 
expected in starburst regions, where there is strong heating of the 
dust by the young stars. 

In our model, the fraction of star formation occurin g in the 
burst mode increases with redshift (see lBaugh et al.ll2005l) . so the 
average IMF with which stars are being formed shifts from being 
close to a solar neighbourhood IMF at the present day to being very 
top-heavy at high redshift. In this model, 30% of star formation 
occured in the burst mode when integrated over the past history of 
the Universe, but only 6% of the current stellar mass was formed 
in bursts, because of the much larger fraction of mass recycled by 
dying stars for the top-heavy IMF. We note that our predictions for 
the IR and sub-mm luminosities of starbursts are not sensitive to 
the precise form of the top-heavy IMF, but simply require a larger 
fraction of m ~ 5 — 2QMq stars relative to a solar neighbourhood 
IMF 

We note that the galaxy formation model in this paper, un- 
like some other recent semi-analytical models, does not include 
AGN feedback. Instead, the role of AGN feedback in reducing the 
amount of gas cooling to form massive galaxies is taken by super- 
winds driven by supernova explosions. These superwinds eject gas 
from galaxy halos, reducing the mass of hot gas and hence also 
the rate of gas cooling in halos. Th e first semi-analytical model to 
include AGN feedback was that of Granato et ai](|2004)), who in- 
troduced a detailed model of feedback from QSO winds during the 
formation phase of supermassive black holes (SMBHs), with the 
aim of explaining the co-evolution of the spheroid al components 
of galaxies and their SMBHs. The predictions of the lGranato et al.l 
model for number counts and redshift distributions in the IR have 
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been computed by Isilva et aL I ( I2OO5I) using the GRASIL spec- 
trophoto metric model, and co mpared to ISO and Spitzer data. How- 
ever, the^ Granato etal.l ( l2004h model has the limitations that it does 
not include the merging of galaxies or of dark halos, and does not 
treat the formation and evolution of galactic disks. More recently, 
several semi-analytical models have been published which propose 
that heating of halo gas by relativistic jets from an AGN in an op- 
tically inconspicuous or "radio" mode can balance radiative cool- 
ing of gas in h igh-mass halos, t hus suppre ssing hot accretion of gas 
onto galaxies I Bower et al.l2006 : Croton e^lj2006l : ICattaneo et al.l 
l2006ljMonaco et al.l2007l : Lagos et aU20d 8l). However, these AGN 
feedback models differ in detail, and all are fairly schematic. None 
of these models has been shown to reproduce both the observed 
number counts and and the observed redshifts of the faint sub-mm 
galaxies. 

The effects of our superwind feedback are qualitatively quite 
similar to those of the radio-mode AGN feedback. Both superwind 
and AGN feedback models contain free parameters, which are ad- 
justed in order to make the model fit the bright end of the ob- 
served present-day galaxy luminosity function at optical and near- 
IR wavelengths. However, since the physical mechanisms are dif- 
ferent, they make different predictions for how the galaxy lumi- 
nosity function should evolve with redshift. Current models for the 
radio-mode AGN feedback are very schematic, but they have the 
advantage over the superwind model that the energetic constraints 
are greatly relaxed, since accretion onto black holes can convert 
mass into energy with a much higher efficiency than can supernova 
explosions. We will investigate the predictions of models with AGN 
feedback for the IR and sub-mm evolution of galaxies in a future 
paper. 

2.2 GRAS I L model for stellar and dust emission 

For each galaxy in our model, we compute the spectral energy dis- 
tribut i on using the spectro photometric model GRASIL jSilva et aU 
ll998l ; lGranato et al.ll200(]|) . GRASIL computes the emission from 
the stellar population, the absorption and emission of radiation by 
dust, and also ra dio emission (therm al and synchrotron) powered 
by massive stars jBressan et al.ll2002h . 

The main features of the GRASIL model are as follows: 

(i) The stars are assumed to have an axisymmetric distribution in a 
disk and a bulge. 

(ii) The SEDs of the stellar populations are calculated separately for 
the disk and the bulge, according to the distribution of stars in age 
and metallicity that is obtained from the corresponding star forma- 
tion and chemical enrich ment histories. We use the population syn- 
thesis model described in lBressan et al.l ( Il998h . which is based on 
the Padova stellar evolution tracks and Kurucz model atmospheres. 
This model is able to reprod uce fairly well the in tegrated UV prop- 
erties of Globular Clusters jChavez et al.ll2009D and the observed 
Spitzer IRS spectra and mid-IR colours of ellipticals in the V irgo 
and Coma clusters jBressan et alj|200^ : IClemens et al.ll2009h . i.e. 
old stellar populations likely at the two extremes of the metallicity 
range of stellar systems. At intermediate and young ages, it com- 
pares well with r ecent observations of LMC star clusters (see e.g. 
iMolla et al.ll2009l) 

(iii) The cold gas and dust in a galaxy are assumed to be in a 2- 
phase medium, consisting of dense gas in giant molecular clouds 
embedded in a lower-density diffuse component. In a quiescent 
galaxy, the dust and gas are assumed to be confined to the disk, 
while for a galaxy undergoing a burst, the dust and gas are con- 
fined to the spheroidal burst component. 



(iv) Stars are assumed to be bom inside molecular clouds, and then 
to leak out into the diffuse medium on a timescale tcsc- As a result, 
the youngest and most massive stars are concentrated in the dustiest 
regions, so they experience larger dust extinctions than older, typ- 
ically lower-mass stars, and dust in the clouds is also much more 
strongly heated than dust in the diffuse medium. 

(v) The extinction of the starlight by dust is computed using a ra- 
diative transfer code; this is used also to compute the intensity of 
the stellar radiation field heating the dust at each point in a galaxy. 

(vi) The dust is modelled as a mixture of graphite and silicate grains 
with a continuous distribution of grain sizes (varying between SA 
and 0.25 ^m), and also Polycyclic Aromatic Hydrocarbon (PAH) 
molecules with a distribution of sizes. The equilibrium temperature 
in the local interstellar radiation field is calculated for each type 
and size of grain, at each point in the galaxy, and this information 
is then used to calculate the emission from each grain. In the case of 
very small grains and PAH molecules, temperature fluctuations are 
important, and the probability distribution of the temperature is cal- 
culated. The detailed spectrum o f the PAH em i ssion is obtained us- 
in g the PAH cross-s ections from lLi & Draind l l200lh . as described 
in lVega et ai]( l2005h . The grain size distribution is chosen to match 
the mean dust extinction curve and emissivity in the local ISM, and 
is not varied, except that the PAH abundance in molec ular clouds 
is ass umed to be 10"'' of that in the diffuse medium jVega et al] 
I2OO5I) . 

(vii) Radio emission from ionized gas in HII regions and syn- 
chrotron radiation from relativistic electrons acceler ated in super- 
nova r emnant shocks are calculated as described in lBressan et al.l 
( I2OO2I) . 

The output from GRASIL is then the complete SED of a 
galaxy from the far-UV to the radio (wavelengths lOOA < A < 
Im). The SED of the dust emission is computed as a sum over 
the different types of grains, having different temperatures depend- 
ing on their size and their position in the galaxy. The dust SED is 
thus intrinsically multi-temperature. GRASIL has been shown to 
give an excellent match to the observed SEDs of galaxies of all 
types, from passi ve systems through to ULIRGs jSilva et al]| 19981 ; 
IVegaetalJl2008l) . We show some example SEDs from the com- 
bined GALFORM + GRASIL code in Fig|T] 



3 EVOLUTION OF THE GALAXY LUMINOSITY 
FUNCTION IN THE FAR-IR 

We begin by showing what the model predicts for the evolution 
of the galaxy luminosity function (LP) at the far-IR wavelengths 
which will be probed by Herschel, before moving on to predic- 
tions for more immediately observable quantities in the next sec- 
tion. Fig. |2] shows how the predicted luminosity function evolves 
from z = to 2 = 6. The top-left panel shows the LF of the 
total mid-l-far-IR luminosity Lir, defined as the integral over the 
galaxy SED from 8 to lOOO/xm. This range includes almost all of 
the dust emission from a galaxy, but hardly any of the stellar emis- 
sion. The total IR LF evolves strongly with redshift, increasing in 
characteristic luminosity and number density from z = to 2: ~ 2, 
having a plateau from 2 ~ 2 to z 4, and then gradually declin- 
ing at higher redshifts. This evolution results from the combined 
effects of the evolution of the galaxy SFR distribution and the shift 
in the dominant mode of star formation from quiescent (with a nor- 
mal IMF) at low redshift to burst (with top-heavy IMF) at higher 
redshift. - this is discussed further in ^ Starting from very high 
redshift, the cosmic star formation density increases to a peak at 
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Figure 2. Predicted evolution of the galaxy LF in the far-IR. The LFs are shown at redshifts z = 0, 1, 2, 3, 4, 6, with different redshifts shown by different 
colours, as indicated by the key. (a) Total IR (8-1000/im) luminosity, (b) Rest-frame 100/tm. (c) Rest-frame 250/tm. (d) Rest-frame 500/tm. 



2 --^ 2 — 3, driven by the build-up of dark matter halos in which 
gas is able to cool and form stars, and then declines down to z = 0, 
driven by the declining efficiency of gas cooling in halos as their 
masses become even larger. Over the redshift range 2; ~ — 3, the 
characteristic luminosity increases by a factor ~ 8, while the num- 
ber density of LIRGs (Luminous Infrared Galaxies, defined to have 
LiR = 10" L0) increases by a factor ~ 10, and the number den- 
sity of ULIRGs (Ultra-luminous Infrared Galaxies, defined to have 
LiR = 10^^ 1/0 ) increases by a factor ~ 300. The other panels in 
Fig. |2] show the predicted evolution of the LF at rest-frame wave- 
lengths of 100, 250 and 500^im (calculated through the PACS and 
SPIRE bandpasses). The evolution of the LF at these wavelengths 
is qualitatively similar to that of the total IR LF. The 100/^m LF 
looks very similar to the total IR LF, both in the form of the evo- 
lution and in normalization (when vL^ is used as the luminosity 
variable), reflecting the fact that the far-IR SEDs peak around this 
wavelength. At longer wavelengths, the characteristic luminosity is 
fainter (in terms of vLi,), reflecting the decline in the SED long- 
wards of the peak; the degree of evolution is also somewhat less (a 
factor 5 in characteristic luminosity at 500/xiti over the redshift 



range 2: O — 3), reflecting a shift in the average far-IR SED shape 
to somewhat warmer dust temperatures at higher redshifts. 

4 NUMBER COUNTS AND REDSHIFT DISTRIBUTIONS 
4.1 Number counts 

We now move on to make predictions for the quantities which can 
be measured most directly in cosmological galaxy surveys with 
Herschel. The planned cosmological surveys will concentrate on 
imaging in the three PACS bands (centred at 70, 100 and 160^m) 
and the three SPIRE bands (centred at 250, 350 and 500^m); Her- 
schel spectroscopy will be limited to fewer and brighter sources. 
The simplest observable quantity is the number counts per solid 
angle of galaxies as a function of flux, Sv, which we plot in the 
fovmdN/d\nS^. 

An observational estimate of the number counts in the SPIRE 
bands has alrea dy been been made using the BLAST b alloon- 
bome telescope jPevlin et al.ll2009l : |Patanchon et alj|2009l) . so we 
first compare our model predictions with those data in Fig. [3] The 
BLAST results are for a sing le 9 deg^ field. The BLAST maps have 
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log(S„/mJy) log(S„/mJy) log(S„/mJy) 

Figur e 3. Predicted galaxy differential number counts in the SPIRE bands (blue curves) compared to observational estimates from BLAST jPatanchon et alj 
12003) . shown in black. The observational estimates were obtained by modelling the P(D) distiibution, rather than identifying individual sources. The points 
show the "best estimates", while the error bars and shaded regions show the 68% confidence ranges. 



two times worse angular resolution than Herschel and are also 
quite noisy, which leads to serious problems with confusion, in- 
completeness and Eddington bias. As discussed by the BLAST au- 
thors, this means that a direct determination of the counts by iden- 
tifying individual sources is not feasible. They therefore estimated 
the counts using a "P(Z))" analysis, fitting a piecewise power-law 
model to the distribution of pixel brightnesses in the maps at each 
wavelength. The results of this are shown in the figure, with the 
points showing the best estimate of the counts at each of the nodes 
of the model fit, and the error bars and shaded region showing 
the 68% confidence region. Note that the BLAST P{D) analysis 
assumed that the sources were not clustered. Including clustering 
could both change the best estimate values, and increase the er- 
ror bars on the estimated counts. We see that the predicted counts 
agree well with the BLAST estimates at faint fluxes, but are a little 
higher at bright fluxes, especially at the shorter wavelengths. How- 
ever, given the aforementioned uncertainties in the BLAST results, 
and the limited sky coverage, these differences should not be taken 
as conclusive at this stage, and a definitive test must await the re- 
sults from Herschel. 

Next, in Fig. |4] we show the predicted differential number 
counts in the PACS and SPIRE imaging bands, for the full range of 
fluxes and wavelengths covered by Herschel, with each panel cor- 
responding to a different wavelength. In each panel, the blue curve 
shows the predicted total counts, while the green and red curves 
show the separate contributions to these from quiescent galaxies 
and ongoing bursts. (A note about our terminology: by "bursts" we 
mean any galaxy with ongoing star formation in the burst mode, 
while by "quiescent" we mean all other galaxies, whether under- 
going star formation only in the disk mode, or completely passive 
with no current star formation. However, completely passive sys- 
tems contribute very little to the luminosity functions and number 
counts at far-IR wavelengths.) We see that quiescent galaxies dom- 
inate the counts at brighter fluxes in all bands, while the bursts tend 
to dominate at faint fluxes, reflecting the increasing domination of 
the starburst mode of star formation at higher redshifts. 

We can use our model to predict the flux levels at which 
sources should become confused in the different Herschel bands. 
Source confusion happens when multiple sources are separated in 
angle by less than the angular resolution of the telescope, and so 
appear merged together in images. Since the number of sources in- 
creases with decreasing flux, confusion sets a lower limit to the flux 



A 


dpWHM 


^conf 


Iconf 


250 


290 




[arcsec] 


[mjy] 








70 


5.2 


0.24 


0.88 


1.15 


2.53 


100 


7.7 


1.6 


0.95 


0.80 


1.76 


160 


12 


9.5 


1.28 


0.72 


1.57 


250 


19 


22 


1.78 


0.70 


1.64 


350 


24 


21 


2.17 


0.93 


1.93 


500 


35 


19 


2.58 


1.21 


1.93 



Table 1. Predicted confusion limits for Herschel imaging. A is the wave- 
length, BpwHM is the angular resolution of PACS or SPIRE imaging at 
that wavelength, S^onf is the predicted flux at the confusion limit, and 
Iconf is the slope of the differential source counts at that flux. 2:50 and zgo 
are the predicted median and 90-percentile redshifts for galaxies brighter 
than Sconf- 



at which one can still identify separate sources in an image (regard- 
less of the integration time used). Confusion will be a serious prob- 
lem in deep cosmological surveys with Herschel due to the rela- 
tively poor angular resolution of the telescope (compared to optical 
telescopes). This confusion limit depends both on the angular reso- 
lution of the telescope and on the actual surface density of sources 
per solid angle as a function of flux. We estim ate the confusion limit 
using the source density criterion ( C ondonI 1 1 9741 ; Fvkisanen et al] 
I2OOII) : if the telescope has a FWHM (full width at half maximum) 
beamwidth of Ofwhm, we define the effective beam solid angle as 
^beam = (tt/ (4 In 2) ) 6l|^y = 1.136l|,^-^jvj, and then define 
the confusion limited flux Sconf to be such that A^(> Sconf) = 
l/(A/'i,eami^6eam), whcrc A'^(> S) is the numbcr per solid angle of 
sources brighter than flux S. This estimate ignores any clustering of 
the sources. We choose Mbcam = 20 for the number of beams per 
source at the confusion limit, w hich g ives similar resul ts to more 
detailed analyses dViiisanen et al.l200ll : lDole et al.il2004h . We have 
calculated confusion limits in the different Herschel bands using 
our predicted number counts, together with the assumed beamsizes 
given in Table[T] The predicted confusion limited fluxes are given in 
the table and plotted as vertical black dashed lines in Fig.|4] Since 
Herschel images will be essentially diffraction-limited, confusion 
sets in at a lower source density at longer wavelengths, which typi- 
cally implies a brighter flux. 

We also indicate in Fig. |4] by vertical and horizontal dashed 
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Figure 4. Galaxy differential number counts in the six Herschel bands. We show three different curves for our standard model - solid blue: total counts 
including dust extinction and emission; solid red: ongoing bursts; solid green: quiescent galaxies. The vertical dashed black lines show the estimated confusion 
limit for the model at each wavelength, calculated as described in j|4. 1| The vertical and horizontal dashed coloured lines show the flux and area limits for 
some planned Key Project surveys, as indicated in the key. 
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Figure 1. Example SEDs from GALFORM+ GRASIL. (a) Quiescently star- 
forming disk galaxy at 2 = 0. (b) Starburst at 2 = seen one e-folding 
time after start of burst. The blue lines show the "stellar" emission (which 
includes emission from dust in AGB star envelopes, and thermal and syn- 
chrotron radio emission from the ISM), both with (solid lines) and without 
(dashed lines) dust extinction. The red lines show emission from interstellar 
dust, both the total (solid line) and the separate contiibutions from molecu- 
lar clouds (dashed line) and the diffuse ISM (dotted line). The shaded green 
region shows the wavelength range 60-600/im covered by the Herschel 
imaging bands. 



coloured lines the regions of flux and surface density that are 
planned to be probed by the main cosmological galaxy surveys 
with Herschel. For each survey, the vertical line indicates the nom- 
inal flux limit Smin set by integration time and signal-to-noise 
(ignoring possible source confusion), and the horizontal line indi- 
cates the minimum surface density of sources that can be probed 
given the solid angle A of the survey, which we estimate as 
{dN/d\n S)min = 1/A. We consider here the followitig four 
planned surveys, which are all Herschel Key Programme^. The 
first three are deep surveys, while the last is a shallower wide-area 
survey. The deep surveys have various tiers, but for simplicity we 



here consider only the deepest blank-field tier in each survey, since 
this sets the limit for how far down in luminosity the survey can 
probe at each redshift, except for HERMES, where we also con- 
sider a shallower tier. 

GOODS-HerscheB 

This survey using PACS and SPIRE will have 2 tiers, "ultradeep" 
and "superdeep", and will include the deepest imaging in the PACS 
100 and 160/im bands of any of the cosmological key programmes. 
We consider here the ultradeep tier, which is in the GOODS-S field. 
PACS Evolutionary Probe fPEPfl' 

This deep imaging survey using PACS, which is coordinated with 
the HERMES survey, will include blank fields covering a range of 4 
in limiting flux and 50 in area, mainly in the 100 and 160/im bands. 
We consider here the deepest blank field, which is the GOODS-S 
field. 

Herschel Multi-tiered Extragalactic Survey (HERMES^ 
This survey, using both PACS and SPIRE, and coordinated with 
PEP, will be the largest of the cosmological surveys in terms of 
observing time. It will include 6 tiers of blank field surveys, cov- 
ering a range of 6 in limiting flux and 400 in area. We consider 
here the deepest tier. Level- 1 (the CDFS/GOODS-S field), and also 
a shallower tier, Level-5 (the XMM, ELAIS-Nl-SCUBA, Bootes- 
SCUBA2, EGS-SCUBA2, CDFS and Lockman fields), which is in- 
termediate in area and in depth in the SPIRE bands between Level- 1 
and the shallower ATLAS survey. 

Herschel ATLaQ 

This survey, using PACS and SPIRE, will be the shallowest of the 
cosmological surveys, but will cover by far the largest area. There 
is only a single tier. 

Table|2]lists the basic parameters (wavelengths, areas and flux 
limits ignoring confusion) for these planned surveys (or tiers within 
surveys) which we will show in subsequent figures. These surveys 
are plotted in Fig. |4] These survey parameters may be modified 
once the inflight performance of the telescope is known. 

It can be seen from Fig. |4] that, taken together, the Her- 
schel surveys at each wavelength (excepting 70/im) should con- 
tain galaxies covering a huge range (~ 10^ — 10^) in flux. At 
70/im, the planned Herschel surveys probe only slightly deeper 
than surveys already carried out using Spitzer. We show in the 
panel the regions of flux and source density probed by two of 



these Spitzer surveys , SWIRq^ ( A = 49deg^, Sr, 



iLonsdale et ai] I2OO3I) and FIDEljl (A _ = 0.5deK^ Srm n 



ISmJy 



3.3mJv?Dickinson & FIDEL team' 20071: iHuvnh et al .''2007'). The 
Herschel PEP GOODS-S survey will probe about two times fainter 
than Spitzer but still well above the predicted Herschel confusion 
limit. At 100/im, the PEP survey is predicted to be at the confusion 
limit, while the GOODS-HERSCHEL survey will be below it. Con- 
fusion is predicted to be worst in the 160/im deep surveys, where 
the sensitivity of the GOODS-HERSCHEL survey will be 10 times 
below the predicted confusion limit, and that for PEP GOODS-S 
6 times below it. At the longer SPIRE wavelengths (250, 350 and 
500/im), the deepest blank field tier (LI) of the HERMES survey 
will be 4-5 times fainter than the predicted confusion limit, while 
L5 tier will be at or slightly below confusion. The ATLAS survey 



|http://herschel.esac.esa.int/Key-ftogrammes.shtml| 



http://herschel.esac.esa.int/Docs/KPOT 

http://www.mpe.mpg.de/ir/Research/PEP 

http://astronomy.sussex.ac.uk/~sjo/Hermes 

http://herschel.esac.esa.int/Docs/KPOT 

http://swire.ipac.caltech.edu/swire 

http://irsa.ipac.caltech.edu/data/SPITZER/FIDEL 
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Table 2. Basic parameters for the surveys (or tiers within surveys) being modelled. The area A is the total value for that tier, and may include multiple 
fields. Limiting fluxes Smin are 5(7 and ignore confusion. For HERMES Level-5 some of the flux limits vary between different fields, so we give the values 
corresponding to the majority of the survey area. Ng^i is the predicted total number of galaxies in that tier, and 250 and 290 are their predicted median and 
90-percentile redshifts. 



is predicted to be safely above the confusion limit (by factors of 
2-AQ) at all wavelengths. 

Various techniques will allow one to probe observation- 
ally sources fainter than the formal confusion limit, the princi- 
pal ones being pixel brightness distributions, gravitational lensing 
and multi-wavelength analysis. In pixel brightness distribution (or 
P{D)) methods, one uses the distribution of pixel brightness to 
make statistical inferences about the slope and amplitude of the 
source counts at and below the confusion limit, without trying to 
identi fy the individual sources responsible (e.g. IPatanchon et al.l 
I2OO9I) .' In gravitational lensing, one uses the gravitational magni- 
fication of the images of faint background galaxies by a foreground 
galaxy cluster to reduce confusion effects. Since both source fluxes 
and areas get multiplied by the magnification A > 1, the cumu- 
lative source counts transform as A^(> S) = {1/ A)N(,{> S/A), 
where No{> S) is the unlensed distribution. If the source counts 
have the power-law form dN/d\n S oc S^'' , then the intrinsic (un- 
lensed) flux at the confusion limit is reduced by a factor A^^^~' 
due to lensing. Table [T] lists the source count slopes predicted at 
the confusion limit for the different Herschel bands. It can be seen 
that the predicted slope 7con/ at the confusion limit increases with 
increasing wavelength, implying that gravitational lensing can po- 
tentially allow one to probe further below the Herschel confusion 
limits at shorter wavelengths. For example, for a gravitational mag- 
nification A — 10, the effective confusion limit is lowered by fac- 
tors of 0.09-0.4 as the wavelength increases from 100 to 500/im. 
Finally, in multi-wavelength analysis, one combines images at dif- 
ferent wavelengths having different angular resolutions. Variants 
of this include multi-wavelength priors, where one starts from a 
source list obtained from higher angular resolution data at some 
other wavelength, and tries to extract fluxes for individual confused 
sources at the target wavelength, and multi-wav elength stacking , 
where one tries to m easure mean fluxes only (e.g. lDole et alj200q : 
iMarsden et alj2009l) . 



4.2 Redshift distributions 

Having identified sources in images and measured their counts as 
a function of flux, the next key step observationally is to measure 
their redshifts (either spectroscopically or using photometric red- 
shift methods), and construct redshift distributions in different flux 
ranges. We therefore consider the model predictions for redshift 
distributions next. 

Firstly, in Fig. |5] we show the predicted median redshift (to- 
gether with the 10% and 90% percentiles) as a function of flux in 
each of the PACS and SPIRE bands. As in Fig. ID we also indi- 
cate the flux at the confusion limit and at the different Key Project 
survey limits by vertical dashed black and coloured lines respec- 
tively. The predicted median and 90% percentile redshifts for each 
survey tier are also given in Table |2l We see from the figure that 
the median redshifts for sources at the confusion limit are around 
Z50 ~ 1 — 1.5, with the highest values at the shortest and longest 
Herschel wavelengths. If one assumes that the confusion limit can 
be completely circumvented, then of the surveys listed, HERMES- 
LI at 500/im will probe to the highest median redshift (2:50 ~ 1.8). 
If one assumes instead that confusion sets a hard limit, then the 
highest median redshift is reduced to 250 ~ 1.4, again achieved in 
the HERMES-Ll survey at 500/im. 

We now look at the redshift distributions in more detail. Fig.|6] 
shows the predicted redshift distributions for galaxies at the confu- 
sion limits listed in Table [T] In this figure, the reshift distributions 
have all been normalized to unit area under the curve to allow easier 
comparison of their shapes. We have included the redshift distribu- 
tion at the 70/im confusion limit for completeness, even though 
none of the planned surveys at 70/im will go this faint. Apart from 
500/im, all of the redshift distributions peak at quite modest red- 
shifts, z ~ 0.4 — 0.8, although there is a tail of objects to z ~ 2. 
The redshift peak gets broader with increasing wavelength, until at 
the longest wavelength, 500/im, it splits into two peaks, with the 
main peak at 2: « 1.4 and a smaller peak at 2: « 0.2. This effect 
at 500/im is a manifestation of the negative k-correction, whereby 
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Figure 5. Median and percentiles of predicted redsliift distributions as functions of flux in the six Herschel bands. The solid blue lines show the median and 
the dashed lines the 10% and 90% percentiles for galaxies at each flux. The vertical dashed black line shows the estimated confusion limit for the model at 
each wavelength. The coloured Unes show the flux limits for some planned Key Project surveys, as indicated in the key. 
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redshifting of the SED combined with the negative slope of the dust 
SED longwards of the peak counteracts the dimming due to the in- 
creasing luminosity distance, and makes higher redshift galaxies 
more easily visible than lower redshift galaxies. This effect is al- 
ready well k nown from longer- w avelength sub-mm observations at 
850/im (e.g. lHughes etai]|l998h . 

The next figure. Fig. [7] shows the predicted redshift distribu- 
tions at different wavelengths for all of the deep survey tiers listed 
in Table[2] In this figure, the redshift distributions are normalized to 
the expected number of galaxies in the survey area, to allow an eas- 
ier estimate of the number of galaxies predicted in different redshift 
ranges for each of the surveys. Fig. [8] examines the redshift distri- 
butions for these deep surveys in more detail, showing the separate 
contributions of quiescent and bursting galaxies to the total redshift 
distributions for selected survey tiers at particular wavelengths. We 
show two survey tiers (GOODS-Herschel and HERMES-Ll) and 
three wavelengths (100, 250 and 500 ^ra) to illustrate the general 
behaviour. In all cases, we see that the quiescent galaxies dominate 
the distribution at low redshifts and the bursts at high redshifts, re- 
flecting the higher luminosities of the bursts. We also see that the 
bursts become more dominant overall at longer wavelengths in the 
deep surveys, due to the effects of the negative k-correction. 

Finally, in Fig.|9l we show the predicted redshift distributions 
for the shallower but wider-area ATLAS survey. Out of the 5 wave- 
lengths in this survey, the largest number of galaxies should be seen 
at 250/im. The median redshift for galaxies brighter than the flux 
limit is also largest at this wavelength (2:50 ~ 0.4), with 20, 000 
galaxies at 2: > 1.3 and ~ 1000 galaxies at z > 2. The redshift dis- 
tribution is broader at the longer wavelengths (250-500 /xm) com- 
pared to the shorter wavelengths (100 and 160/im), again as a result 
of the negative k-correction at longer wavelengths. 



5 PHYSICAL PROPERTIES OF Herschel GALAXIES 

Having presented predictions for directly observable quantities 
(fluxes and redshifts) in the previous section, we now move on to 
the predicted physical properties of the galaxies detected in Her- 
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Figure 8. Redshift distributions for galaxies brighter than the flux limits of 
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Figure 7. Redshift distributions in planned deep blank-field surveys, (a) GOODS-Herschel. (b) PEP GOODS-S. (c) HERMES LI (CDFS). (d) HERMES L5. 
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schel surveys as a function of flux and redshift. We focus here on 
four properties of central physical importance: the total IR lumi- 
nosity, LiR, the star formation rate, SFR, the stellar mass, Mstar, 
and the host dark halo mass, Mhaio- 

In Fig. [To] we plot the median total IR luminosity Ljr (inte- 
grated over the wavelength range 8-lOOO^m) against flux, for three 
different Herschel bands (one in each panel), for galaxies at one of 
six different redshifts (2 = 0.25, 0.5, 1, 2, 3, 4, as indicated by the 
different colours shown in the key). We have chosen wavelengths 
of 100, 250 and 500/xni to be representative of the six imaging 
bands. The "error bars" on each line show the 10-90% range of the 
distribution at each flux and redshift. We have also plotted vertical 
lines showing the nominal flux limits of the different Key Project 
surveys discussed in the previous section, as well as the predicted 
confusion limit at each wavelength. As mentioned previously, Lir 
essentially measures the total luminosity of dust emission from a 
galaxy. For galaxies with significant recent star formation, Lir is 
powered mostly by far-UV radiation from massive young stars, and 



it thus provides a tracer of the dust-obscured star formation rate for 
high mass (m > 5A/0) stars. The actual conversion factor between 
Lir and the total SFR (integrated over all stellar masses) depends 
both on the fraction of the far-UV light from young stars which is 
absorbed by dust (which is typically high) and on the IMF. We see 
from Fig. [To] that at each Herschel wavelength and each redshift, 
there is an approximately linear relationship between Lir and the 
flux in that band with only modest scatter (~ 0.2dex). This sim- 
ply reflects the fact that the Herschel bands directly probe the rest- 
frame far-IR wavelengths which dominate Lir, and that the shape 
of the far-IR SED shows only modest galaxy-to-galaxy variation at 
a given redshift, with only a weak dependence on Lir. The zero- 
point of this relation between Lir and flux obviously depends on 
the SED shape and on the effects of the luminosity distance and 
the k-correction; the zero-point changes with redshift less at longer 
wavelengths (over the range z = 0.5 — 4, it increases by 2.5 dex 
at lOO^m and by ~ 1.5 dex at 500/im), reflecting the effect of the 
negative k-correction at the longer wavelengths. 
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Figure 9. Redshift distributions at 100, 160, 250, 350 and SOO^tm at the flux 
limits of the ATLAS survey. Different wavelengths are shown by different 
colours, as indicated in the key. The distributions are all normalized so that 
the integral over the distribution is equal to the expected number of galaxies 
in the survey area. Note that in this case the log of the redshift distribution 
is plotted. 



One of the main goals of Herschel will be to measure the evo- 
lution of the cosmic density of dust-obscured star formation, and 
for this purpose it is interesting to know which Herschel wave- 
length will probe to the faintest Lir at each redshift, since the lat- 
ter will determine what fraction of the total IR luminosity density 
is actually resolved into identified objects at that redshift. We see 
from Fig.llOlthat if source confusion sets a hard limit to identifying 
individual sources, then surveys at lOO/im will probe to the faintest 
Lir, ranging from ~ lO^°-^/i"^L0 at z = 0.5 to ~ IQ^'^-^h'^LQ 
at z = 4. On the other hand, if all of the planned Key Project sur- 
veys manage to resolve objects down to their nominal flux limits 
(even where these are below confusion), then the GOODS-Herschel 
survey at 160/im will probe faintest (down to Lir ^ l(f '^h~^LQ 
at 2 = 0.5 and ~ lO^^^/i"^!/© at z = 4). Of the surveys in the 
SPIRE bands at 250-500/.im, those at will probe down to 

the lowest Lir whether confusion can be circumvented or not, ex- 
cept at the highest redshifts, z > 2 — 3, for which the 500^m band 
becomes more sensitive due to the negative k-correction effect. 

Fig- El is similar to Fig.fTO] except that the SFR rather than 
Lir is plotted against flux. As already described, the relation be- 
tween Lir and SFR depends on the fraction of UV light from 
young stars absorbed by dust and on the IMF. In particular, in 
our model, Lir/ SFR is about 4 times larger for star formation 
in bursts with the top-heavy {x = 0) IMF compared t o the stars 
forming quiescently in disks with the iKennicuttI ( Il983h IMF. The 
proportion of star formation associated with the burst mode on av- 
erage increases with increasing Lir, and for this reason the relation 
between SFR and flux is shallower than a linear proportionality. 
The scatter is also somewhat larger than for the L/h-Aux relation. 
We see from the figure that, at the confusion limit, Herschel sur- 
veys should probe down to SFRs ~ Ih^^ Mc-yyv^^ z — 0.5 and 
~ 10^ MQyr~^ at z — 4, the lowest limits being achieved at 
100//m. The planned Key Project surveys may improve on this by 
factors ~ 2 if they can get below the confusion limit. 

Fig.[l2]is analogous to the previous two figures, but now with 
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Figure 10. Median total IR (8-lOOO/^m) luminosity, Lir, vs flux at 100, 
250 and 500 fim for galaxies selected at different redshifts, shown by dif- 
ferent colours as indicated in the key. The "error bars" on the lines show the 
10-90% range. The vertical dashed black line shows the predicted confusion 
limit at each wavelength, while the vertical dashed coloured lines show the 
nominal flux limits for different planned surveys, as indicated in the key. 
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Figure 11. Median total SFR vs flux at 100, 250 and 500 /im for galaxies 
selected at different redshifts. The different lines are as described for Fig llOl 



Figure 12. Median stellar mass v.v flux at 100, 250 and 500 /^m for galaxies 
selected at different redshifts. The different lines are as described for Fig llOl 
For clarity, we have introduced small horizontal offsets between the lines 
plotted for different redshifts. 
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the stellar mass Mstar plotted against flux. In this case, the relation 
is far from linear, and has much larger scatter than for Lm or SFR. 
This is not surprising, since the Herschel flux is proportional to the 
emission from dust heated mostly by massive young stars, while 
the stellar mass includes stars of all ages and masses. Since the 
galaxies in our model have complex star formation histories, there 
is no simple relation between the current star formation rate and 
the total mass of stars formed over the history of the galaxy. The 
galaxies found in Herschel surveys should have stellar masses over 
the range Mstar ~ 10^ — A/q, with a weaker dependence 

on redshift at a given flux than is typically seen for the SFR. 

Finally, in Fig. [13] we show the relation between the host dark 
halo mass Mhaio and flux in the Herschel bands. Here we see that 
the median halo mass depends much more weakly on flux or red- 
shift (and with a larger scatter) than either Mstar or SFR, especially 
at the longer wavelengths. This reflects the fact that in our model 
the relation between far-IR luminosity and halo mass is even more 
indirect than for the stellar mass, especially due to the dominance 
of transient bursts at the higher luminosities. This produces the rel- 
atively flat trend of median halo mass with Herschel flux. The weak 
dependence on redshift is because halos of a given mass on aver- 
age host higher SFRs at higher z, which compensates for the larger 
luminosity distance. We see from Fig. [T3] that the galaxies found 
in the Key Project cosmological surveys should typically have halo 
masses M^aio ~ 10^^h~^ Mq. This will have important implica- 
tions when we consider the clustering of Herschel galaxies in ijs] 

We have also investigated the dependence of the stellar bulge- 
to-total mass ratio B/T and the cold gas mass Mgas on flux in the 
Herschel bands. For brevity, we only show results for the 250/im 
band in Fig. [H We find that the 10-90% range for B/T covers 
nearly the whole possible range ^ B /T ^ 1 at most fluxes and 
redshifts of interest in the planned Herschel cosmological surveys, 
i.e. there is no strong preference for one morphological type over 
another The exception to this is at low fluxes and low redshifts, 
where most sources have B/T < 0.5. This reflects the fact that 
the galaxies found in these surveys should be a mixture of quies- 
cently star-forming disk galaxies and starbursts triggered by galaxy 
mergers, and even though the most luminous galaxies in the far-IR 
are predicted to be bursts, these can be triggered by either major 
or minor mergers, producing a bursting galaxy which can be ei- 
ther bulge or disk dominated. There is a trend for the median B/T 
to increase with flux at a given redshift, presumably reflecting an 
increase in the fraction of bursts triggered by major mergers. We 
illustrate these trends in the top panel of Fig. [14] For the cold gas 
mass Mgas, shown in the lower panel of Fig. [14] we find a corre- 
lation with flux which is weaker than linear, and also fairly weakly 
dependent on redshift. We predict that the galaxies found in the 
Herschel cosmological surveys should typically have gas masses 
~ lO^''ft~^Af0, which implies that many of them should have CO 
emission from their molecular gas which is de tectable by current 
telescopes (e.g. lSolomon & Vanden Bou3l2005l) . 
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Obtaining multi-wavelength data complementary to that from Her- 
schel itself will be crucial for achieving the science goals of the cos- 
mological surveys. Such data will be essential for obtaining iden- 
tifications and accurate positions of the sources, for obtaining ac- 
curate redshifts, and for learning more about the physical nature of 
these galaxies. It is therefore of interest to see what the model pre- 
dicts for the observability of Herschel sources at other wavelengths. 



Figure 13. Median halo mass \'.v flux at 100, 250 and 500 fim for galaxies 
selected at different redshifts. The different lines are as described for Fig llOl 
For clarity, we have introduced small horizontal offsets between the lines 
plotted for different redshifts. 
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Figure 15. Median fluxes or magnitudes at other wavelengths for galaxies selected at 100/jm. (a) GALEX NUV. (b) r-band. (c) IRAC 3.6^m. (d) MIPS 24/jm. 
(e) SCUBA 850/im. (f) 1.4GHZ. The triangles next to the the y-axes in these panels indicate the magnitude or flux limits for the different surveys which are 
discussed in the text. For clarity, we have introduced small horizontal offsets between the lines plotted for different redshifts. 
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Figure 16. Median fluxes or magnitudes at other wavelengths for galaxies selected at 250/im. (a) GALEX NUV. (b) r-band. (c) IRAC 3.6^m. (d) MIPS 24/jm. 
(e) SCUBA 850/im. (f) 1.4GHZ. The triangles indicate the flux limits for different surveys, as in Fig. [15] For clarity, we have introduced small horizontal 
offsets between the lines plotted for different redshifts. 
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Figure 17. Median fluxes or magnitudes at other wavelengths for galaxies selected at 500/jm. (a) GALEX NUV. (b) r-band. (c) IRAC 3.6/^m. (d) MIPS 24/jm. 
(e) SCUBA 850/im. (f) 1.4GHZ. The triangles indicate the flux limits for different surveys, as in Fig. [15] For clarity, we have introduced small horizontal 
offsets between the lines plotted for different redshifts. 



20 Lacey et al. 



—\ 1 1 r 




J I I L_ 



1 2 
log(S,(250)/mJy) 




8.5 



z=2 

confusion z = 
ATfJiS z=4 
HDRMES L5 
HEiRMES LI 



_I I I I I L 



_I I I L_ 



1 2 
log(S,(250)/mJy) 



Figure 14. Median bulge-to-total stellar mass ratio B/T and cold gas mass 
vs flux at 250 /tm for galaxies selected at different redshifts. The differ- 
ent lines are as described for Fig llOl For clarity, we have introduced small 
horizontal offsets between the lines plotted for different redshifts. 



In Figs. [15] and [T7] we show predicted fluxes or AB magni- 
tudes at wavelengths from the UV to the radio, plotted against 
the Herschel flux at either 100 or 250 or 500^m. Specifically, 
we consider the GALEX NUV filter (centred at A = 0.23/xm), 
the SDSS r-band (centred at A = 0.62/im), the Spitzer 3.6 and 
24/im bands, the SCUBA 850Atm band, and the VLA 1.4GHz 
band. The observer-frame NUV, 24^m, 850/im and 1.4GHz bands 
in the main trace recent star formation, while the observer-frame 
3.6/im band traces older stars, and the observer-frame r-band traces 
older stars at low redshift but younger stars at high redshift (when 
it corresponds to the rest-frame UV). The predicted median fluxes 
at 0.23, 24 and 850/im and at 1.4 GHz generally track the fluxes 
in the Herschel bands fairly well, but with zero-points which can 
depend strongly on redshift, depending on the wavelength. This is 
particularly the case for the GALEX NUV band, for which the flux 
drops rapidly with redshift at 2: > 1.5 when the band falls short- 
ward of the Lyman break in the rest frame. (Note that we include 
the effects of absorption by the intervening IGM when we calculate 
the UV and optical magnitudes.) The redshift dependence of the ze- 



ropoints is generally least at 1 .4 GHz and 24^m. As expected, the 
correlation with the Herschel flux is much weaker for the r and 
3.6/im bands, which are more sensitive to stellar mass than to the 
SFR. 

For comparison, we note the approximate flux or magnitude 
limits for some of the main surveys at these wavelengths (the flux 
limits for these surveys are also indicated by triangles along the 
y-axis in Figs. ll5ll7t : 

GALEX NUV: All-sky Imaging Survey (AIS) ruAB = 20.8; 
Medium Imaging Survey (MIS) rriAB = 22.7; Deep Imaging Sur- 
vey (DIS) TUAB = 24.4 i Morrissev et alj2007h . 

22.2 jAbazaiian et"al 



SDSS r: SDSS Legacy Survey tuab 



2009); Subaru Deep Field rriAs = 27.8 (Rc) l lKashikawa et al 
2004); Hubble Ultra -Deep Field uiab = 30.1 (V or i) 
( lBeckwithet"ai]|2006h ; Pan-STARRS PSl tuab = 24.1 in the 
Stt survey and rriAB = 26.9 in the Medium Deep Survey (MDS) 

ItChambers 2006). 

Spitzer 3.6^im: SWI RE S„ = 3.7MJy jUonsdale et al.ll2004l) ; 
GOODS = 0.6/iJy (IDickinson et al.ll2003l) . Also the WISE 
satellite w ill survey the whole sky to Sv = 120/xJy at a wavelength 

of 3.3Atm jWright et alj2007l) . 

Spitzer 24fj,m: SWI RE = lOO/xJy jLonsdale et aLll2004l) ; 
GOODS = 70^iJy dChary et alj[2004t) . Also WISE will sur- 
vey the whole sky to Sv = 2600/iJy at a wavelength of 23/xm 

dWright et al.ll2007l) . 

SCUBA 850fj.m: HDF 5"^ = 2mJy (Hughes etal]|l998h ; 
SHADES 5"^ = 8mJy jMortier et alj |2005) SASSjQ 5^. = 

150mJy. 

1.4 GHz: N VSS Sy = 2.5mJy jCondon et al.|[T99g) ; FIRST 
Sv = ImJv teeckeretalj I1995I) ; Phoenix Sy = 60/jJy 
llHopkinsetai]|2003l) ; SSA 13 Sv = 20/iJy jpomalont et"al] 
I2OO6I) . 



7 UNVEILING THE COSMIC STAR FORMATION 
HISTORY 

Since the discovery of the far-IR background by COBE, it has been 
known that the bulk of star formation over the history of the Uni- 
verse has been obscured by dust. One of the primary goals of Her- 
schel is to resolve the far-IR background into individual sources 
and hence determine the amount of dust-obscured star formation 
at different cosmic epochs. How well Herschel will be able to do 
this depends on the distribution of the total IR emissivity, e/jj, (i.e. 
the mean luminosity density per comoving volume) over sources of 
different luminosities and over redshift, and how far down in total 
IR luminosity, Lm, Herschel surveys are able to probe at different 
redshifts. This section is devoted to investigating the implications 
of our models for this key issue. 

We start by showing in the top left panel of Fig.[T8]the cosmic 
star formation history predicted by our model, both the total star 
formation density and the separate contributions to this from ongo- 
ing starbursts and from quiescently star-forming galactic disks (this 
was earlier shown in Baugh et al. 2005). The total SFR density in- 
creases by a factor ~ 6 from z — to 2 — 2.5, and then declines 
very gradually to higher z. We see that the quiescent mode of star 
formation dominates the SFR density at redshifts z < 3 probed by 
Herschel, while the burst mode dominates at higher redshifts. 



http ://w w w.j ach.hawaii . edu/JCMT/ surveys/sassy/ 
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Figure 18. (a) Cosmic SFR history, sliowing tlie total SFR density (blue) and the separate contributions from bursts (red) and quiescent disks (green). The 
solid lines show the total SFR density integrated over all stellar masses, while the dashed lines show the SFR density in massive stars (m > SMq). The dotted 
lines show the rate of build-up of stellar mass in long-lived stars and remnants, after allowing for recycling of gas to the ISM. (b) The time integral of SFR 
density from t = up to t{z). The line colours and styles are identical to (a), (c) Evolution of luminosity density in rest-frame UV (with and without dust, 
shown by the solid and dashed lines) and in the mid/far-IR. (d) The same as in (c), except plotted on a linear scale against cosmic time. 



We also show by dashed hnes in the same panel the star for- 
mation density in massive stars only, which we define as stars 
with masses m > 5Mq, which have lifetimes < 1 x lO^yr. 
We choose this mass range because, when integrated over the 
whole galaxy population, such stars dominate the UV emissivity 
from galaxies, and also dominate the heating of the dust which 
powers most of the mid- and far-IR emission from galaxies. The 
two IMFs in our model (assumed to cover the stellar mass range 
0.15 < m < I2OM0) have very different fractions of their initial 
stellar mass in high mass stars: we find /(m > 5Mq) = 0.24 
for the Kennicutt IMF assumed for quiescent star formation, and 
/(m > 5i\f0) — 0.96 for the top-heavy a: = IMF assumed for 
bursts. (For comparison, /(m > 5Mq) — 0.22 for a Salpeter IMF 
covering the same mass range.) The SFR density for massive stars 
therefore evolves more strongly than that for all stars, increasing 
by a factor ~ 15 from 2; = to a peak at 2; ~ 3. For massive star 
formation, the burst mode already dominates the quiescent mode at 
z > 1. Finally, the dotted lines in the top panel show the rate of 



build-up of stellar mass in long-lived stars and stellar remnants, af- 
ter accounting for the mass returned to the ISM by dying stars. The 
fraction of the initial stellar mass returned in this way is called the 
recycled fraction, and depends on the IMF, having values of 0.41 
for the Kennicutt IMF and 0.91 for the a; = IMF. We see that, 
after allowing for recycling, the quiescent mode of star formation 
dominates the build-up of stellar mass even at the highest redshifts 
plotted. 

We show in the top right panel of Fig. [TS] the time integral 
of the SFR density from t = up to redshift z for the different 
components. Focussing on the integral up to 2; = 0, we see that, in- 
tegrated over the history of the Universe, 31% of all star formation 
is predicted to occur in the burst mode, but at the present day this 
fraction is only 4.8%. For high mass (m > 5A/0) stars, 64% of 
star formation happened in the burst mode over the history of the 
universe, while the present-day fraction is 17%. The burst mode is 
thus more important for high mass star formation. However, after 
accounting for recycling of mass from dying stars, we find that the 
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fraction of the present-day stellar mass produced by the burst mode 
is only 6.4%. 

In the lower panels of Fig. [TS] we show how the model SFR 
density evolution translates into the evolution of the cosmic emis- 
sivity from galaxies in the UV (defined here as the integral over 
the wavelength range 0.01 — 0.3/im) and the total mid/far-IR (in- 
tegrated over the range 8 — lOOO/^m, as previously). The UV emis- 
sivity is plotted both with and without dust extinction. The unex- 
tincted UV and total IR emissivities both increase by factors ~ 10 
from 2: = to z = 3, and then remain approximately constant up 
to z = 6. This difference in redshift dependence compared to that 
for the total SFR density is because both the UV and total IR emis- 
sivities are powered mostly by massive stars (m > SM©), and the 
burst mode of star formation produces a larger fraction of such stars 
than the quiescent mode. We see that the effect of dust extinction 
on the UV emissivity is predicted to be very large - integrated over 
the history of the Universe, 90% of the UV energy is predicted to 
be absorbed by dust. This fraction increases from 77% at z = to 
87% at 2 = 1 and 92% at z — 6. This emphasizes how essential 
it is to measure the cosmic evolution of the total IR emissivity in 
order to measure directly the cosmic SFR history, free of uncertain 
observational estimates of UV dust extinction. The total IR emis- 
sivity is seen to be very similar to the unextincted far-UV emis- 
sivity at all redshifts, which just reflects the fact that most of the 
far-UV radiation is absorbed by dust, and that heating of the dust 
by longer wavelength (optical and near-IR) radiation from stars is 
a minor contribution to the total when integrated over the whole 
galaxy population. (The fraction of the IR emissivity due to heat- 
ing by longer wavelength radiation is predicted to be around 20% 
integrated over the history of the Universe, dropping from 40% at 
2: = to only 20% at 2 = 1.5 and 10% at 2 = 3.) 

We next consider the fraction of the total IR emissivity, em, 
which is produced by galaxies brighter than total luminosity Lir 
for different redshifts. This is obtained by integrating over the total 
IR luminosity function, which was shown in Fig. [2] The results are 
shown in the top and bottom panels of Fig. [19] The top panel shows 
the fraction of the IR emissivity in galaxies brighter than Lir as 
a function of Lir for different redshifts, while the bottom panel 
shows the same fraction as a function of redshift for a number of 
different minimum luminosities. We see from the top panel that the 
value of Lir above which 50% of the IR emissivity is contributed 
increases with increasing 2, from 3 x lO^'^h'^ Lq at z — to 
2 X 10^^ Lq at 2 = 3, due to the general brightening in the 
luminosity function with 2 discussed in ij3l This evolutionary effect 
partly compensates for the effect of increasing luminosity distance, 
which generally means that the minimum Lir for which we can 
detect galaxies in surveys increases with 2. 

The middle panel in Fig.[T9]shows the fraction of star forma- 
tion occuring in galaxies brighter than Lir, as a function of Lir 
for different redshifts. We show this both for the total SFR inte- 
grated over all stellar masses (solid lines) and for the high-mass 
(m > 5Mq) SFR (dashed lines). As a consequence of the bright 
end of the total IR LF being dominated by bursts with a top-heavy 
IMF, galaxies brighter than a given Lir account for a larger frac- 
tion of the high-mass SFR than of the total SFR. For example, at 
2 = 3, galaxies with Lir > 2 x W^^ Lq account for 50% of 
the high-mass SFR but only 36% of the total SFR. 

The critical question is: what fraction of the total IR emissiv- 
ity will different Herschel surveys be able to resolve at different 
redshifts? To answer this, we plot in the top panel of Fig. 1201 the 
fraction of the total IR emissivity eiR produced by galaxies brighter 
than some flux Sv as a function of redshift. Each of the six Herschel 
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Figure 19. (a) Fraction of total IR luminosity density from galaxies brighter 
than Lir at different redshifts. (b) Fractions of total SFR (solid lines) and 
high-mass SFR (dashed lines) from galaxies brighter than Lir at different 
redshifts. (c) Fraction of total IR luminosity density from galaxies brighter 
than Lir as function of redshift. 
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Figure 20. (a) Fraction of IR luminosity density from galaxies brighter than 
expected Herschel survey flux limits at different wavelengths. The soHd 
lines correspond to the faintest nominal flux limit at each wavelength (ig- 
noring confusion) from Table |2] while the dashed lines correspond to the 
confusion limit at each wavelength from Table [T] (b) Similar to (a), but 
showing the fraction of the time integral of the IR luminosity density up to 
redshift z which is in galaxies brighter than flux Sy in a Herschel band. 



imaging bands is plotted in a different colour, and for each band we 
show the result for two different flux limits: the solid line shows the 
result for the faintest planned flux limit at that wavelength ignor- 
ing confusion, taken from the list of surveys in Table (2] while the 
dashed line shows the result at the estimated confusion limit, taken 
from Table[T] The lower panel of Fig.|20]instead shows the fraction 
of the time-integrated emissivity Uir{z) = Jq*'^' e/ij dt which is 
resolved at a given flux limit. We see from this figure that if source 
confusion is ignored (or can be circumvented), then the GOODS- 
Herschel Ultra-deep survey at 160^m (with a planned flux limit of 
0.9mJy) should resolve the largest fraction of the total IR emis- 
sivity into sources at all redshifts. Integrated over all redshifts, a 
survey at this flux limit would resolve 49% of all IR emission into 
sources. On the other hand, if confusion sets a hard limit, then the 
GOODS-Herschel Ultra-deep survey at 100/^m (with a confusion 
limit of 1.6mJy) should instead resolve the largest fraction of the 
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Figure 21. (a) Mean clustering bias fe as a function of Ljjj for galaxies at 
different redshifts, as indicated in the key. (b) Clustering length ro obtained 
from the bias and power spectrum as function of L jij for different redshifts. 

IR emissivity at all redshifts, and a total fraction of 27% when in- 
tegrated over all redshifts. (In principle, a confusion-limited sur- 
vey at would resolve a larger fraction, but no such survey 
is planned). Note that the question we have asked (and answered) 
here is different from asking what fraction of the present-day cos- 
mic IR background (GIB) is resolved into sources at different fluxes 
and wavelengths. The energy density in the GIB at 2: = is given 
by the integral fg° em/il + z) dt, which differs from the time- 
integrated emissivity Uir{z — 0) by a factor 1/(1 + z) in the 
integrand due to redshifting of the photon energies. 

In summary, the planned surveys with Herschel should be able 
to resolve into sources around 30-50% of the total IR dust emission, 
and thus a similar fraction of the dust-obscured massive star forma- 
tion, over the history of the Universe. 



8 CLUSTERING OF GALAXIES IN THE FAR-IR 

The final topic we consider is the predicted clustering of galax- 
ies detected in the far-IR by Herschel. In the GALFORM model. 
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the clustering of galaxies is determined (on large scales) by that 
of their host dark matter halos, and (on small scales) by the spa- 
tial distribution of galaxies within dark matter halos. Observational 
measurements of clustering for galaxies of different classes and lu- 
minosities can thus provide robust information about the masses of 
the dark halos hosting these galaxies. Such information is difficult 
to obtain in other ways, particularly for high redshift galaxies (for 
example, measurements of galaxy circular velocities do not probe 
the radii where most of the halo mass is, and so only weakly con- 
strain halo masses). Clustering measurements can thus provide an 
essential test of the link between mass (mostly in dark matter) and 
light, which provides a fundamental test of the galaxy formation 
model. 

We limit ourselves here to a few predictions for large-scale 
clustering, and defer a more detailed analysis of clustering proper- 
ties to a future paper. On large scales, where the density fluctuations 
for the dark matter are still in the linear regime, the clustering of 
galaxies and halos can be described by a linear bias factor b, such 
that^(r) — fo^^dm(?'), where ^(r) is the 2-point corelation function 
for the galaxies or halos, and (,dmij') is that for the dark matter. We 
can calculate the bias, bhaio{M, z), for halos as a func tion of halo 
mass a nd redshift using the analytical approximation of lSheth et al.l 
( l200lh . We can then calculate the bias, bgai{Lii{, z), for galaxies 
of a specified IR luminosity as a mean of the halo bias, weighted 
by the distribution of host halo masses for galaxies of this lumi- 
nosity and redshift. The latter information is provided by the semi- 
analytical model. We show the results of this calculation in the top 
panel of Fig. |2T] We see from this that the clustering bias is pre- 
dicted to depend only weakly on total IR luminosity over the range 
10'^ < LiR < 10^'^-^ h-'^LQ, but to increase steeply at the highest 
luminosities. The weak dependence on IR luminosity stems from 
the weak correlation of IR luminosity with host halo mass, which 
was seen in Fig. [13] The bias at a given IR luminosity is seen to 
increase gradually with redshift, from 6 1 at z = 0, to 6 ~ 2 at 
2 = 3, for the luminosity range where b is roughly constant. 

By combining our analytical estimate of the bias with an 
estimate of the correlation function, (,dm{r,z), of the dark mat- 
ter, we can calculate the correlation function of the galaxies, 
^gai{r,LiH,z) = 6g„,(L/B,2:)^dm(r,2). This estimate win be 
valid on scales large enough that the linear bias approximation 
is valid. We calculate the dark mat ter correlation fun ction using 
the approximate analytical model of lSmith et al.l | |2003|) . which al- 
lows for non-linear effects in the dark matter evolution. From this, 
we can calculate the correlation length, ro, which we define by 
^(ro) = 1, and which provides a simple, directly measureable, 
characterization of how strongly a particular population of galax- 
ies is clustered. This approach to calculating ro is valid provided 
& ^ 1, since then (,dm{ro) ^ 1, and the dark matter is still ap- 
proximately in the linear regime on the scale ro. We show our 
predictions for ro in the lower panel of Fig. [2T] We see from 
this that, as for the bias, ro is almost independent of luminos- 
ity over the range 10^ < Ljr < W^^ '^h'^LQ, but increases 
steeply for Lir > 10"'^/i-^Lq. At a given luminosity, ro de- 
creases with increasing redshift, but typically only gradually. For 
10^ < LiR < lO"-^fe~^L0, ro decreases from ro ~ 5/i~^Mpc 
at z = to ro ~ 3/i~^ Mpc at 2 = 3. 

IVieroetai] ( l2009l) have estimated the clustering of galaxies in 
the SPIRE bands from the BLAST observations by measuring the 
angular power spectra of the total intensity maps. From this, they 
estimate a bias 6 ~ 4 for the galaxies responsible for the far-IR 
background, larger than the values typical in our model. However, 
their estimate of the bias relies on a model for the source redshift 



distribution, and assumes that the bias is independent of both red- 
shift and IR luminosity. The latter assumption seems unrealistic in 
the light of our own predictions for the bias shown in Fig. l2U a). 
We therefore do not make a more detailed comparison with their 
results. 

Finally, in Fig. (22] we show simulated images of a slice 
through the universe at 2 = 1, 100^~^Mpc wide by 10/i~^Mpc 
thick. These have been obtained by combining our s emi-analytical 
mode l with the Millennium dark matt er simulatio n ( Springel et al.l 
using the same method as in lOrsi et alj ( |2008|) . In these 
images, the dark matter distribution is shown in green, while the 
positions and luminosities of the galaxies are shown as coloured 
blobs. For making these images, we have calculated the far-IR lu- 
minosities o f galaxies using the si mplified dust emission model 
described in [Gonzalez et al.l ( |2010() . since running the GRASIL 
dust code separately on each galaxy in the Millennium simulation 
would not have been feasible computationally. In a future paper, 
we will present images calculated using a more ac curate approxi- 
matio n to GRAS I L based on Artificial Neural Nets jAlmeida et al.l 
I2OO9I) . The top two panels show galaxies selected based on their 
total IR luminosities, brighter than lO"/i~^L0 and lO^^/i'^L© 
respectively in the left and right panels. The lower left panel shows 
galaxies selected to have 100/im fluxes brighter than 2mJy, which 
is similar to the planned flux limit in the PEP GOODS-S sur- 
vey. Finally, the lower right panel shows galaxies selected accord- 
ing to their dust-extincted rest-frame far-UV luminosities, with 
Mas(1500A) ~ 5 log ft < -21.2. This absolute magnitude limit 
has been chosen because it corresponds to about the same SFR in 
a completely unextincted galaxy as does Ljr = 10^^ Lq in a 
completely extincted galaxy. Comparing the top left and lower right 
panels illustrates how incomplete surveys for star-forming galaxies 
can be if they use only observations in the rest-frame far-UV and 
ignore the far-IR. 



9 SUMMARY AND CONCLUSIONS 

We have used a detailed model of hierarchical galaxy formation 
and of the reprocessing of starlight by dust to make predictions 
for the evolution of the galaxy population at the far-infrared wave- 
lengths (60-670/im) which will be probed by observations with 
the Hersche! Space Observatory. We calculated galaxy formation 
in the framework of the ACDM model using the GALFORM semi- 
analytical model, which includes physical treatments of the hier- 
archical assembly of dark matter halos, shock-heating and cool- 
ing of gas, star formation, feedback from supernova explosions and 
photo-ionization of the IGM, galaxy mergers and chemical enrich- 
ment. We computed the IR luminosities and SEDs of galaxies using 
the GRASIL multi- wavelength spectrophotometric model, which 
computes the luminosities of the stellar populations in galaxies, and 
then the reprocessing of this radiation by dust, including radiative 
transfer through a two-phase dust medium, and a self-consistent 
calculation of the distribution of grain temperatures in each galaxy 
based on a local balance between heating and cooling. 

Our galaxy formation model incorporates two different IMFs: 
quiescent star formation in galaxy disks occurs with a normal solar 
neighbourhood IMF, but star formation in starbursts triggered by 
galaxy me rgers happens with a top-heavy x — IMF. In a previ- 
ous paper teaugh et alj|2005h . we found that the top-heavy IMF in 
bursts was required in order that the model reproduces the observed 
number counts of the faint sub-mm galaxies detected at 850 fim, 
which are typically ultra-luminous starbursts at 2 --^ 2, with total IR 



Predictions for Herschel 25 




Figure 22. Images of a simulated slice of the universe IQOh ^Mpcwideand \0h ^Mpcthickat^ = 1. Each panel shows the same slice, with the dark matter 
density plotted in green, and with the galaxies plotted as coloured blobs, the blob size increasing with luminosity or flux. Each panel shows galaxies selected at 
a different wavelength and/or lumninosity/flux. (a) L/^ > lO" /i^^L©. (b) L/^ > lO^^h'^ Lq. (c) S„(100Aim) > 2mjy. (d) Myv - 5ioc//i < -21.2. 



luminosities Lm 10^ — 10 Lq. We subsequently found that 
the same model also reproduces the evolution of the galaxy pop- 
ulation at mid-infrared wavelengths found by Spitzer (Lacev et al. 
1 200 8 ). We have used the same model, with identical parameters, to 
make predictions for Herschel in the present paper. 

We began by showing the predicted evolution of the 
galaxy luminosity function at far-IR wavelengths. This brightens by 
a factor ~ 10 going back from z = to 2: ~ 3, reflecting both the 
evolution of star formation rates in galaxies and the increasing im- 
portance of the top-heavy burst mode with increasing redshift. We 
next presented predictions for galaxy number counts as func- 
tions of flux in the Herschel PACS and SPIRE imaging bands (cov- 
ering the wavelength range 70-500/im). We calculated the confu- 
sion limits for the Herschel bands, and found that source confusion 
is likely to be a serious problem for all of the deepest cosmological 
surveys planned (PEP GOODS-Herschel and HERMES) (except 
at lOfim). The number of faint galaxies which can be resolved in 
these surveys will depend dramatically on whether or not the con- 
fusion limit can be circumvented, e.g. by using multi-wavelength 
data. We also investigated the predicted redshift distributions in 



these deep cosmological surveys and in the wide-area ATLAS sur- 
vey. We found that the deep surveys should reach median redshifts 
~ 1 — 1.8, depending on the wavelength and on whether it is possi- 
ble in practice to probe sources fainter than the confusion limit. The 
highest median redshift (1.4-1.8, depending on confusion) should 
be attained in the HERMES survey at 500^m. For the ATLAS sur- 
vey, the median redshift should be 0.2-0.4, with the highest value 
at 250/im. At the faintest fluxes and highest redshifts, the galaxy 
source population is predicted to be dominated by starbursts. 

Following on from the predictions for source counts and red- 
shifts, in ^we showed what the model predicts for some of the 
basic physical properties of galaxies detected in different Her- 
schel bands as functions of flux and redshift. As expected, at each 
redshift there are nearly linear corelations between the fluxes in 
different bands and the total IR luminosity Lm (integrated over 
the wavelength range 8-lOOO^m). On the other hand, the relation 
between the total star formation rate SFR and flux shows more 
non-linearity, principally due to the different IMFs assumed in 
starburst galaxies (which dominate at high luminosity) and qui- 
escent galaxies (which dominate at low luminosity). The deep- 
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est surveys should resolve galaxies with Lir > 10^^ Lq at 
z — I and Lin 10^^ Lq at z = 2, corresponding to SFRs 
> 10 — lOO/i^^ Afoyr"^. The stellar and dark halo masses of Her- 
schel galaxies show much weaker correlations with flux and red- 
shift, and with much more scatter, in large part because of the role 
of starbursts triggered by galaxy mergers. In the deep surveys at 
2 = 1 — 2, the typical //eric/ieZ-detected galaxy should have a stel- 
lar mass ~ 10"'/i"^MQ and a halo mass ~ W^h~^MQ, and a 
typical gas mass ~ lQ^^h~^ Mq. Finally, the morphologies of the 
galaxies detected in the Herschel cosmological surveys should be 
quite mixed, rannging from highly bulge-dominated to highly disk 
dominated systems, reflecting our assumption that starbursts can be 
triggered by both major and minor galaxy mergers. 

Since our model is a multi-wavelength model, we used it in 
^to predict what should be the fluxes at other wavelengths (rang- 
ing from the far-UV to the radio) of galaxies detected in Herschel 
surveys. Follow-up data at other wavelengths will be essential both 
to determine redshifts of Herschel sources and to investigate their 
physical properties. Specifically, we presented predictions for the 
GALEX NUV band, the SDSS r band, the Spitzer 3.6 and 24^m 
bands, the SCUBA band and finally the 1.4GHz radio flux. 

One of the primary goals of Herschel will be to unveil the 
history of dust-obscured star formation in the Universe. Therefore, 
in ^ we presented the predictions of our model for the cosmic 
star formation history and for the evolution of the cosmic emissiv- 
ity of galaxies in the UV and the mid/far-IR. Our model predicts 
that, over the history of the Universe, about 90% of the UV radi- 
ation from massive young stars has been reprocessed by dust into 
the mid/far-IR wavelength range. We used our model to investigate 
what fraction of the total energy emitted by dust heated by stars 
over the history of the Universe should be resolved into galaxies by 
different planned Herschel surveys. We find that the fraction of the 
total IR emission resolved in this way should be ~ 50% if individ- 
ual sources can be resolved all the way down to the nominal survey 
flux limits set by signal-to-noise, but only ~ 30% if source con- 
fusion sets a hard limit. Of the currently planned surveys, those at 
100 or 160/im should be able to resolve the largest fraction of the 
time-integrated IR emission. This then implies that Herschel should 
resolve a similar fraction of the massive star formation over the his- 
tory of the Universe (roughly m > SM© , since these stars are re- 
sponsible for most of the UV heating of dust grains). Lower-mass 
stars make only a small contribution to powering the mid/far-IR 
emission from dust, so determining the total SFRs from Herschel 
observations relies on extrapolating to lower masses (m < SMq) 
based on an assumed IMF. If, as assumed in our model, the IMF is 
different in different types of galaxy (quiescent vi- starburst), then 
estimating total SFRs from IR data becomes much more compli- 
cated and uncertain, but estimates of the SFRs in high-mass stars 
should be much more robust. 

Finally, in ijs] we briefly investigated the predicted clustering 
of Herschel galaxies. We found that the typical galaxies in Her- 
schel cosmological surveys should have a modest clustering bias 
(fe ~ 1 — 2) relative to the dark matter, with correlation lengths 
ro ~ 3 — 5/!.~^Mpc (in comoving units), except at the highest lu- 
minosities {Lir > 10^^ Lq). Measurements of the clustering 
in Herschel surveys will provide an essential test of the relation be- 
tween galactic star formation rates and halo masses predicted by 
our galaxy formation model. 
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